function [trackResults, channel]= tracking_BladeRF(fidi,fidq, channel, settings)
% Performs code and carrier tracking for all channels.
%[trackResults, channel] = tracking(fid, channel, settings)
%   Inputs:
%       fid             - file identifier of the signal record.
%       channel         - PRN, carrier frequencies and code phases of all
%                       satellites to be tracked (prepared by preRum.m from
%                       acquisition results).
%       settings        - receiver settings.
%   Outputs:
%       trackResults    - tracking results (structure array). Contains
%                       in-phase prompt outputs and absolute starting 
%                       positions of spreading codes, together with other
%                       observation data from the tracking loops. All are
%                       saved every millisecond.
%--------------------------------------------------------------------------
%                           SoftGNSS v3.0
% 
% Copyright (C) Dennis M. Akos
% Written by Darius Plausinaitis and Dennis M. Akos
% Based on code by DMAkos Oct-1999
%
% Mod via Xiaoyeyimier
% Mod @ 2021-10-11
% Crop : CAUC
% Adjust to fit BladeRF to work
%--------------------------------------------------------------------------
%% Initialize result structure ============================================
%---------- Channel status %-----------------------------------------------
trackResults.status         = '-';                                         % No tracked signal, or lost lock
trackResults.absoluteSample = zeros(1, settings.msToProcess);              % 1 ms per code 
trackResults.codeFreq       = inf(1, settings.msToProcess);                % Frequency of CA code
trackResults.carrFreq       = inf(1, settings.msToProcess);                % Frequency of the tracked carrier wave:
%---------- Outputs from the correlators (In-phase) -----------------------
trackResults.I_P            = zeros(1, settings.msToProcess);              % The branch I prompt
trackResults.I_E            = zeros(1, settings.msToProcess);              % The branch I early
trackResults.I_L            = zeros(1, settings.msToProcess);              % The branch I late
%---------- Outputs from the correlators (Quadrature-phase) ---------------
trackResults.Q_E            = zeros(1, settings.msToProcess);              % The branch Q prompt
trackResults.Q_P            = zeros(1, settings.msToProcess);              % The branch Q early
trackResults.Q_L            = zeros(1, settings.msToProcess);              % The branch Q late
%---------- Loop discriminators -------------------------------------------
trackResults.dllDiscr       = inf(1, settings.msToProcess);
trackResults.dllDiscrFilt   = inf(1, settings.msToProcess);
trackResults.pllDiscr       = inf(1, settings.msToProcess);
trackResults.pllDiscrFilt   = inf(1, settings.msToProcess);
%---------- Copy initial settings for all channels ------------------------
trackResults = repmat(trackResults, 1, settings.numberOfChannels);         % Maximum number of signal processing channels supported
%% Initialize tracking variables  
codePeriods = settings.msToProcess;                                        % For GPS 1 C/A code is 1 ms 

%---------- DLL variables -------------------------------------------------
% Define early-late offset (in chips)
earlyLateSpc = settings.dllCorrelatorSpacing;                              % Ahead or behind by half a code chip
PDIcode = 1e-3;                                                            % Summation interval = 1 [ms]
[tau1code, tau2code] = calcLoopCoef(settings.dllNoiseBandwidth,...         % Calculate filter coefficient values
                                    settings.dllDampingRatio,1.0);         % settings.dllNoiseBandwidth = 2;[Hz] 
                                                                           % settings.dllDampingRatio = 0.7;
                                                                           
%---------- PLL variables -------------------------------------------------
PDIcarr = 1e-3;                                                            % Summation interval = 1 [ms]
[tau1carr, tau2carr] = calcLoopCoef(settings.pllNoiseBandwidth, ...        % Calculate filter coefficient values
                                    settings.pllDampingRatio,0.25);        % settings.pllNoiseBandwidth　= 25;[Hz]
                                                                           % settings.pllDampingRatio = 0.7;
hwb = waitbar(0,'Tracking...');                                            % create a waitbar  
acq_Num = 0;                                                               % channel counter 
for channelNr = 1:settings.numberOfChannels                                % check how many channels acquired
    if (channel(channelNr).PRN ~= 0)                                       % if a satellite is captured
        acq_Num = acq_Num + 1;                                             % counter plus one 
    end
end
debugfig = figure();                                                       % gen a new figure
axCount = 0;                                                               % axes counter 
codeFreq      = settings.codeFreqBasis;                                    % settings.codeFreqBasis = 1.023e6; [Hz]

%% Start processing channels ==============================================
for channelNr = 1:settings.numberOfChannels
    
    % Only process if PRN is non zero (acquisition was successful)
    if (channel(channelNr).PRN ~= 0)
        axCount = axCount + 1;
        % Save additional information - each channel's tracked PRN
        trackResults(channelNr).PRN     = channel(channelNr).PRN;   %%%将捕获的卫星号直接作为跟踪的卫星
        
        % Move the starting point of processing. Can be used to start the
        % signal processing at any point in the data record (e.g. for long
        % records). In addition skip through that data file to start at the
        % appropriate sample (corresponding to code phase). Assumes sample
        % type is schar (or 1 byte per sample) 
        as = fseek(fidi,settings.skipNumberOfBytes + ...
            settings.BytesPerData*(channel(channelNr).codePhase-1),'bof');
        as = fseek(fidq,settings.skipNumberOfBytes + ...
            settings.BytesPerData*(channel(channelNr).codePhase-1),'bof');
        % Get a vector with the C/A code sampled 1x/chip
        caCode = generateCAcode(channel(channelNr).PRN);
        % Then make it possible to do early and late versions
        caCode = [caCode(1023) caCode caCode(1)];

        %--- Perform various initializations ------------------------------

        % define initial code frequency basis of NCO
        
        
        % define residual code phase (in chips)
        remCodePhase  = 0.0;% The initial code shift is 0
        % define carrier frequency which is used over whole tracking period
        carrFreq      = channel(channelNr).acquiredFreq;% The initial carrier frequency is calculated by the capture module
        carrFreqBasis = channel(channelNr).acquiredFreq;
        
        % define residual carrier phase
        remCarrPhase  = 0.0;

        %code tracking loop parameters
        oldCodeNco   = 0.0;
        oldCodeError = 0.0;

        %carrier/Costas loop parameters
        oldCarrNco   = 0.0;
        oldCarrError = 0.0;
        %=== Process the number of specified code periods =================
        for loopCnt =  1:codePeriods
            
%% GUI update -------------------------------------------------------------
            % The GUI is updated every 50ms. This way Matlab GUI is still
            % responsive enough. At the same time Matlab is not occupied
            % all the time with GUI task.
            if (rem(loopCnt, 50) == 0)
                try
                    waitbar(loopCnt/codePeriods, ...
                            hwb, ...
                            ['Tracking: Ch ', int2str(channelNr), ...
                            ' of ', int2str(settings.numberOfChannels), ...
                            '; PRN#', int2str(channel(channelNr).PRN), ...
                            '; Completed ',int2str(loopCnt), ...
                            ' of ', int2str(codePeriods), ' msec']);                       
                catch
                    % The progress bar was closed. It is used as a signal
                    % to stop, "cancel" processing. Exit.
                    disp('Progress bar closed, exiting...');
                    return
                end
            end

%% Read next block of data ------------------------------------------------            
            % Find the size of a "block" or code period in whole samples
            % Update the phasestep based on code freq (variable) and sampling frequency (fixed)
            codePhaseStep = codeFreq / settings.samplingFreq;
            blksize = ceil((settings.codeLength-remCodePhase) / codePhaseStep);% settings.codeLength = 1023;
            % Read in the appropriate number of samples to process this
            % interation 
            [rawSignal_Idata, samplesRead_Idata] = ...
                fread(fidi, blksize, settings.dataType);% settings.dataType should be double here;
            [rawSignal_Qdata, samplesRead_Qdata] = ...
                fread(fidq, blksize, settings.dataType);% settings.dataType should be double here;
           %% for bupt if data
            rawSignal = rawSignal_Idata+ 1j*rawSignal_Qdata;
            rawSignal = rawSignal.';
            %% pulse blanking
            if settings.IsBlanking
                Ex = mean(abs(rawSignal));                                 % Compute the expectation of the Rayleigh distribution
                th = sqrt(2/pi)*sqrt(-2*log(settings.Pa))*Ex;              % Calculate thresholds from expectations
                rawSignal(abs(rawSignal)>th)=0;                            % pulse blanking
            end
            if (samplesRead_Idata ~= blksize)||(samplesRead_Qdata ~= blksize)
                disp('Not able to read the specified number of samples  for tracking, exiting!')
                fclose(fidi);
                fclose(fidq);
                return
            end
%% Set up all the code phase tracking information -------------------------
            % Define index into early code vector
            tcode       = (remCodePhase-earlyLateSpc) : ...
                          codePhaseStep : ...
                          ((blksize-1)*codePhaseStep+remCodePhase-earlyLateSpc);
            tcode2      = ceil(tcode) + 1;
            earlyCode   = caCode(tcode2);
            % Define index into late code vector
            tcode       = (remCodePhase+earlyLateSpc) : ...
                          codePhaseStep : ...
                          ((blksize-1)*codePhaseStep+remCodePhase+earlyLateSpc);
            tcode2      = ceil(tcode) + 1;
            lateCode    = caCode(tcode2);
            % Define index into prompt code vector
            tcode       = remCodePhase : ...
                          codePhaseStep : ...
                          ((blksize-1)*codePhaseStep+remCodePhase);
            tcode2      = ceil(tcode) + 1;
            promptCode  = caCode(tcode2);
            
            remCodePhase = (tcode(blksize) + codePhaseStep) - 1023.0;

%% Generate the carrier frequency to mix the signal to baseband -----------
            time    = (0:blksize) ./ settings.samplingFreq;
            % Get the argument to sin/cos functions
            trigarg = ((-carrFreq * 2.0 * pi) .* time) + remCarrPhase;
            remCarrPhase = rem(trigarg(blksize+1),(2*pi));
            
%% Generate the six standard accumulated values ---------------------------
            % First mix to baseband
            Lo = exp(1j*trigarg(1:blksize));
            BasebandSignal = rawSignal.*Lo;
            iBasebandSignal = real(BasebandSignal); 
            qBasebandSignal = imag(BasebandSignal);

            % Now get early, late, and prompt values for each
            I_E = sum(earlyCode  .* iBasebandSignal);
            Q_E = sum(earlyCode  .* qBasebandSignal);
            I_P = sum(promptCode .* iBasebandSignal);
            Q_P = sum(promptCode .* qBasebandSignal);
            I_L = sum(lateCode   .* iBasebandSignal);
            Q_L = sum(lateCode   .* qBasebandSignal);
%             Ip(loopCnt) = I_P;
%             Qp(loopCnt) = Q_P;
%% Find PLL error and update carrier NCO ----------------------------------
            % Implement carrier loop discriminator (phase detector)
            carrError = atan(Q_P / I_P) / (2.0 * pi);% Phase discrimination factor
            % Implement carrier loop filter and generate NCO command
            carrNco = oldCarrNco+ (tau2carr/tau1carr) * ...
                (carrError - oldCarrError) + carrError * (PDIcarr/tau1carr);
            oldCarrNco   = carrNco;
            oldCarrError = carrError;
            % Modify carrier freq based on NCO command
            carrFreq = carrFreqBasis + carrNco;
            trackResults(channelNr).carrFreq(loopCnt) = carrFreq;
%% Find DLL error and update code NCO -------------------------------------
            codeError = (sqrt(I_E * I_E + Q_E * Q_E) - sqrt(I_L * I_L + Q_L * Q_L)) / ...
                (sqrt(I_E * I_E + Q_E * Q_E) + sqrt(I_L * I_L + Q_L * Q_L));
            % Implement code loop filter and generate NCO command
            codeNco = oldCodeNco + (tau2code/tau1code) * ...
                (codeError - oldCodeError) + codeError * (PDIcode/tau1code);
            oldCodeNco   = codeNco;
            oldCodeError = codeError;
            
            % Modify code freq based on NCO command
            codeFreq = settings.codeFreqBasis - codeNco;
            
            trackResults(channelNr).codeFreq(loopCnt) = codeFreq;

%% Record various measures to show in postprocessing ----------------------
            % Record sample number (based on 8bit samples)
            trackResults(channelNr).absoluteSample(loopCnt) = ftell(fidi);

            trackResults(channelNr).dllDiscr(loopCnt)       = codeError;
            trackResults(channelNr).dllDiscrFilt(loopCnt)   = codeNco;
            trackResults(channelNr).pllDiscr(loopCnt)       = carrError;
            trackResults(channelNr).pllDiscrFilt(loopCnt)   = carrNco;

            trackResults(channelNr).I_E(loopCnt) = I_E;
            trackResults(channelNr).I_P(loopCnt) = I_P;
            trackResults(channelNr).I_L(loopCnt) = I_L;
            trackResults(channelNr).Q_E(loopCnt) = Q_E;
            trackResults(channelNr).Q_P(loopCnt) = Q_P;
            trackResults(channelNr).Q_L(loopCnt) = Q_L;
        end % for loopCnt

        % If we got so far, this means that the tracking was successful
        % Now we only copy status, but it can be update by a lock detector
        % if implemented
% for debug propsose
        axTemp = subplot(acq_Num,2,2*axCount-1);
        plot(trackResults(channelNr).I_P,'b')
        namestr = sprintf('I branch prompt PRN:%2d',trackResults(channelNr).PRN);
        title(namestr);
        axTemp = subplot(acq_Num,2,2*axCount);
        plot(trackResults(channelNr).Q_P,'r')
        namestr = sprintf('Q branch prompt PRN:%2d',trackResults(channelNr).PRN);
        title('Q branch prompt');
        trackResults(channelNr).status  = channel(channelNr).status;        
    end % if a PRN is assigned
end % for channelNr 
debugfig.Position = [100,100,1500,800];
% Close the waitbar
close(hwb)
